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Abstract 

Energy density and energy flux was introduced along Takesue's method. Parti- 
cle energies were localized at particle positions using Dirac delta function. The en- 
ergy density was connected with the energy flux by continuity equation. New method 
was proposed to deal with multi-body interaction. Consequently, the energy current 
between particles could be calculated even if the multi-body interaction appeared. 
Moreover, An application to the molecular dynamics simulation of hydrogen adsorp- 
tion on a graphene with modified Brenner reactive empirical bond order potential was 
demonstrated. 



§1. Introduction 

In classical particle system, the energy current between particles was considered by Lepri, 
Livi and Politi with the approximation of low-A; limit. Takesue produced the energy cur- 
rent between particles from the energy density formed by the Dirac delta function and the 
continuity equation. 2 ' This method is close to the work by Irving and KirkwoodP These 
methods can treat only two-body interaction. However, molecular dynamics (MD) which 
simulates atoms and molecules generally adopt multi-body interaction to represent molec- 
ular structures and chemical properties. In the present paper, the Takesue's method is 
introduced and our new method enables to derive the energy current between particles from 
the multi-body interaction. 

§2. Energy Field by Dirac Delta Function 

We derive the energy density and the energy flux according to Takesue's workP Hamil- 
tonian in many-particle system is given by 

i 

where 7*j(t), Pi(t) and rrii are the z th particle position, momentum and mass, respectively. 
Here, the z-th particle energy ej(t) is defined by 

e i (t) = ^£+u i (t), (2-2) 

where Pi(t) 2 /Irrii and Ui(t) are the z-th particle kinetic energy and interaction energy. The 
particle interaction energies ttj(t) satisfy U(ri(t), r 2 (t), ■ ■ • ) = J2i u ii^)- Using the Dirac 
delta function, the particle energies ei(t) are localized at the particle positions rj(t). Then 
the energy density e(x,t) is defined by 

e(aj,t) = 5^e i (*) < J(aj-ri(t)), (2-3) 

i 

where x = (x, y, z) is space coordinates. 

Because total energy is a conservative quantity, the energy density constructs continuity 
equation 

+ V • jOM) = 0, (2-4) 

with the vector field j(x,t) which is regarded as the energy flux. The energy flux j(x,t) 
should consists of the following two parts. As the z-th particle moves at velocity f,, energy 



2 



current is generated at the i-th particle position r». In addition, interaction transports 
energy between particles by the interaction. From these points of view, the energy flux 
j(x, t) is given by 

j(x, t) = J2 ei{t)ri{t)5{x - n{t)) + ^ ji^ k (t)Xik(x, t), (2-5) 

i i,k>i 

where ji^k(= —jk^i) is the magnitude of the energy current from the i—th particle to the 
fc-th one and Xik{x,t) is the vector field. The first and second terms of the right-hand in 
Eq. (I2 5|) are the energy currents due to the movement of particles and due to interaction 
between particles, respectively. To satisfy Eq. (12-41) . the relation 

V ■ Xl k(x, t) = 5(x - n(t)) - 8{x - r h (f)), (2-6) 

is imposed on the vector field Xik(x,t). Thereby, we obtain 

m k+i 

Thus, if the time derivative of the i—th particle energy de,/dt consists of the elements of 
summation Ylk^v we can re g ar d the elements as the energy current from i—th particle to 
the k— th one jj_>fc(t). To derive ji^k(f) from Eq. (I2-7|) . we must consider the i—th particle 
interaction energy Ui(t). Two-body interaction potential energy is simply allocated to the 
particle energies Ui(f) into equal halves. However, multi-body interaction potential energy 
cannot be divided into the the particle energies Ui(t). 

Determination the particle interaction energy Ui(f) is unnecessary as long as the time 
derivative of the particle interaction energy dui/dt is given. If Hamiltonian system conserves 
total momentum, multi-body interaction potential has to be a function of relative position 
vectors U(f 12(f), ■ ■ ■ ,fik(t),- ■ • , rN-i,N(f)), where fih = ri — rk is the relative position vector 
from the fc-th particle to the i-th one. We note that the multi-body interaction potential 
is not always the function only of the distance between particles = |r^|. 

According to the fact that the multi-body interaction potential is the function of the 
relative position vectors r^, the total differential of the total interaction potential is 

dU(f 12(f), ■ ■ ■ ,f ik (f), ■ ■ ■ ,f N -i )N (f)) = 22 ' & f ik(f) 

i,k>i 

= 5E^-M«)dt. (2.8) 

From this, the time derivative of the particle interaction energy is given by 

dm 1 dU , . , s . . s . , .„ 

if = 2 •(»■•(«) -'»('))«"■ < M ) 

k^i lK 
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Thereby, the time derivative of the i— th particle energy is written by the summation Ylkj^i 
as follows; 

dt dt \2mJ dt 2 i^\ mi m J ' d ^k' 
As a result, the energy current between particles j^ k becomes 

1 ( Pi . Pk\ dU 

2 \rrii m k J dr ik 

The partial derivative dU/dfi k (— —dU/df k i) is the force which acts on the fc-th particle 
due to the variation of the relative position vector in the multi-body interaction. This 
method is effective in arbitrary multi-body potentials. 

§3. Application to Hydrogen Adsorption 

We demonstrate the application of the energy current between particles to the MD sim- 
ulation which dealt with the adsorption of a hydrogen atom on a graphene. The incident 
energy of the hydrogen atom was set to 3 eV. Chemical interaction was represented by 
modified Brenner reactive empirical bond order potentiaP® 



i,j>i 



^(^)-^({r},{^},{e DH })^(r,) 



(3-1) 



where the functions and V^, represent repulsion and attraction, respectively. The func- 



[ij] iv] 

tion bij generates multi-body force where the bond angle 8j ik and the dihedral angle are 
written by cos#? fc = ■ r^/r^r^ and cos 9™ = (r ik x r^) • (r^ x r tj ) /r^r^r^, respec- 
tively. 

The energy current from the hydrogen atom to the graphene j'h = Yli^i 3i-*i an d total 
transport energy = J * j^dt were measured where the index of the hydrogen atom was 1. 
Figure [1] shows j'h and with time. It was understood that the energy of about 2 eV was 
transferred from the hydrogen atom to the graphene as the hydrogen atom is adsorbed. 

§4. Summary 

The energy density and the energy flux was introduced along the Takesue's method. The 
particle energies ej(t) were localized at the particle positions r$(t) using Dirac delta function 
as Eq. (12-31) . The energy density was connected with the energy flux by the continuity 
equation f)2-4p . We proposed the new method to deal with the multi-body interaction. The 
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Fig. 1. (a) the energy current and (b) the total transport energy from the hydrogen atom to the 
graphene. 

total differential of the multi-body interaction potential was expanded by not the particle 
positions Vi but the relative position vectors f^. Consequently, the energy current between 
particles jj-^t) was given by Eq. (12- lip in the multi-body interaction. This method was 
applied to the MD simulation of the hydrogen adsorption on the graphene with modified 
Brenner REBO potential. 
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